/****************************************************************
 * Fulguro
 * Copyright (C) 2004 Christophe Clienti
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 *
 * This library is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
 * Boston, MA 02111-1307, USA.
 ***************************************************************/


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <flgrCoreCopy.h>
#include <flgrCoreMalloc.h>
#include <flgrCoreDispatch.h>
#include "flgrLinearFFT.h"

/*!
 * \defgroup group_fulguro_linear Linear Processing
 * \defgroup group_fulguro_linear_fft Fast Fourier Transform
 * \ingroup group_fulguro_linear
 * Fast Fourier Transform based on fftw3
 * @{
 */


FILE *flgr1d_fft_create_fftw_wisdom_file(int size_x, int DorInv,
					 FLGR_Type type, char *transformtype) {
  char fname[1024];

  FLGR_DECLARE_FUNCTION;

  sprintf(fname,".fftw_wisdom_%dd_%d_%d_%s_%s",1,size_x,DorInv,
	  flgr_get_type_string(type),transformtype);

  return fopen(fname,"wb");
}

FILE *flgr1d_fft_look_for_fftw_wisdom_file(int size_x, int DorInv,
					   FLGR_Type type, char *transformtype) {
  char fname[1024];

  FLGR_DECLARE_FUNCTION;

  sprintf(fname,".fftw_wisdom_%dd_%d_%d_%s_%s",1,size_x,DorInv,
	  flgr_get_type_string(type),transformtype);

  return fopen(fname,"rb");
}



FILE *flgr2d_fft_create_fftw_wisdom_file(int size_x, int size_y, int DorInv,
					 FLGR_Type type, char *transformtype) {
  char fname[1024];

  FLGR_DECLARE_FUNCTION;

  sprintf(fname,".fftw_wisdom_%dd_%dx%d_%d_%s_%s",2,size_x,size_y,DorInv,
	  flgr_get_type_string(type),transformtype);

  return fopen(fname,"wb");
}

FILE *flgr2d_fft_look_for_fftw_wisdom_file(int size_x, int size_y, int DorInv,
					   FLGR_Type type, char *transformtype) {
  char fname[1024];

  FLGR_DECLARE_FUNCTION;

  sprintf(fname,".fftw_wisdom_%dd_%dx%d_%d_%s_%s",2,size_x,size_y,DorInv,
	  flgr_get_type_string(type),transformtype);

  return fopen(fname,"rb");
}







/**************************************** 1D Functions ***************************/











///////////////////////////////////////////////////////////////////
// Macro for different type of flgr1d_fft_forward_cplx
///////////////////////////////////////////////////////////////////
#define FLGR_FORWARD_COMPLEX_FFT_1D(dtype,dtypename,prefix)		\
  FILE *wisdom_file=NULL;						\
  void *bufsrc,*buffft;							\
  prefix##_plan p_forw;							\
  int nx=dat->length;							\
  									\
  FLGR_DECLARE_FUNCTION;						\
									\
  bufsrc = fftw_malloc(nx*sizeof(prefix##_complex));			\
  buffft = fftw_malloc(nx*sizeof(prefix##_complex));			\
									\
  flgr1d_export_raw_##dtype(bufsrc,dat);				\
									\
  wisdom_file=flgr1d_fft_look_for_fftw_wisdom_file(nx,1,dtypename,	\
						   "forward_c2c");	\
  if(wisdom_file!=NULL) prefix##_import_wisdom_from_file(wisdom_file);	\
									\
  p_forw = prefix##_plan_dft_1d(nx, bufsrc, buffft,			\
				FFTW_FORWARD,FFTW_MEASURE);		\
  									\
  if(wisdom_file==NULL) {						\
    wisdom_file=flgr1d_fft_create_fftw_wisdom_file(nx,1,dtypename,	\
						   "forward_c2c");	\
    prefix##_export_wisdom_to_file(wisdom_file);			\
    POST_INFO("fftw wisdom created\n");					\
  }									\
									\
  prefix##_execute(p_forw);						\
									\
  prefix##_destroy_plan(p_forw);					\
									\
  fclose(wisdom_file);							\
  									\
  flgr1d_import_raw_##dtype(datfft,buffft);				\
									\
  fftw_free(bufsrc);							\
  fftw_free(buffft)


void flgr1d_fft_forward_cplx_fgBIT(FLGR_Data1D *datfft, FLGR_Data1D *dat) {
  flgr_no_define_type_function(FLGR_BIT);
}
void flgr1d_fft_forward_cplx_fgUINT8(FLGR_Data1D *datfft, FLGR_Data1D *dat) {
  flgr_no_define_type_function(FLGR_UINT8);
}
void flgr1d_fft_forward_cplx_fgUINT16(FLGR_Data1D *datfft, FLGR_Data1D *dat) {
  flgr_no_define_type_function(FLGR_UINT16);
}
void flgr1d_fft_forward_cplx_fgUINT32(FLGR_Data1D *datfft, FLGR_Data1D *dat) {
  flgr_no_define_type_function(FLGR_UINT32);
}
void flgr1d_fft_forward_cplx_fgINT8(FLGR_Data1D *datfft, FLGR_Data1D *dat) {
  flgr_no_define_type_function(FLGR_INT8);
}
void flgr1d_fft_forward_cplx_fgINT16(FLGR_Data1D *datfft, FLGR_Data1D *dat) {
  flgr_no_define_type_function(FLGR_INT16);
}
void flgr1d_fft_forward_cplx_fgINT32(FLGR_Data1D *datfft, FLGR_Data1D *dat) {
  flgr_no_define_type_function(FLGR_INT32);
}
void flgr1d_fft_forward_cplx_fgFLOAT32(FLGR_Data1D *datfft, FLGR_Data1D *dat) {
  FLGR_FORWARD_COMPLEX_FFT_1D(fgFLOAT32,FLGR_FLOAT32,fftwf);
}
void flgr1d_fft_forward_cplx_fgFLOAT64(FLGR_Data1D *datfft, FLGR_Data1D *dat) {
  FLGR_FORWARD_COMPLEX_FFT_1D(fgFLOAT64,FLGR_FLOAT64,fftw);
}

////////////////////////////////////////////////////////////////////////////
/*! Compute a Complex fft
 *  @param datfft : pointer to FLGR_Data1D complex (spp=2)
 *  @param datsrc : pointer to FLGR_Data1D complex (spp=2)
 *  @returns FLGR_RET_OK, ...
 */
//////////////////////////////////////////////////////////////////
FLGR_Ret flgr1d_fft_forward_cplx(FLGR_Data1D *datfft, FLGR_Data1D *datsrc) {
  FLGR_Ret ret;

  FLGR_DECLARE_FUNCTION;

  if((datsrc == NULL) || (datfft == NULL)) {
    POST_ERROR("Null objects!\n");
    return FLGR_RET_NULL_OBJECT;
  }

  if(datsrc->spp!=2)  {
    POST_ERROR("Only complex images are supported\n");
    return FLGR_RET_NULL_OBJECT;
  }

  if((ret=flgr1d_is_data_same_attributes(datsrc,datfft,__FUNCTION__)) != FLGR_RET_OK) return ret;

  FLGR_DISPATCH_PROCEDURE(datsrc->type, flgr1d_fft_forward_cplx, datfft,datsrc);
}




///////////////////////////////////////////////////////////////////
// Macro for different type of flgr1d_fft_backward_cplx
///////////////////////////////////////////////////////////////////
#define FLGR_BACKWARD_COMPLEX_FFT_1D(dtype,dtypename,prefix)		\
  FILE *wisdom_file=NULL;						\
  void *bufdest, *buffft;						\
  prefix##_plan p_back;							\
  int nx=datfft->length;						\
									\
  FLGR_DECLARE_FUNCTION;						\
									\
  bufdest = fftw_malloc(nx*sizeof(prefix##_complex));			\
  buffft = fftw_malloc(nx*sizeof(prefix##_complex));			\
									\
  flgr1d_export_raw_##dtype(buffft,datfft);				\
									\
  wisdom_file=flgr1d_fft_look_for_fftw_wisdom_file(nx,0,dtypename,	\
						   "backward_c2c");	\
  if(wisdom_file!=NULL) prefix##_import_wisdom_from_file(wisdom_file);	\
									\
  p_back = prefix##_plan_dft_1d(nx, buffft, bufdest,			\
				FFTW_BACKWARD,FFTW_MEASURE);		\
  									\
  if(wisdom_file==NULL) {						\
    wisdom_file=flgr1d_fft_create_fftw_wisdom_file(nx,0,dtypename,	\
						   "backward_c2c");	\
    prefix##_export_wisdom_to_file(wisdom_file);			\
    POST_INFO("fftw wisdom created\n");					\
  }									\
  									\
  prefix##_execute(p_back);						\
  									\
  prefix##_destroy_plan(p_back);					\
  									\
  fclose(wisdom_file);							\
									\
  flgr1d_import_raw_##dtype(datdest,bufdest);				\
									\
  fftw_free(bufdest);							\
  fftw_free(buffft);							\
									\
  return


void flgr1d_fft_backward_cplx_fgBIT(FLGR_Data1D *datdest, FLGR_Data1D *datfft) {
  flgr_no_define_type_function(FLGR_BIT);
}
void flgr1d_fft_backward_cplx_fgUINT8(FLGR_Data1D *datdest, FLGR_Data1D *datfft) {
  flgr_no_define_type_function(FLGR_UINT8);
}
void flgr1d_fft_backward_cplx_fgUINT16(FLGR_Data1D *datdest, FLGR_Data1D *datfft) {
  flgr_no_define_type_function(FLGR_UINT16);
}
void flgr1d_fft_backward_cplx_fgUINT32(FLGR_Data1D *datdest, FLGR_Data1D *datfft) {
  flgr_no_define_type_function(FLGR_UINT32);
}
void flgr1d_fft_backward_cplx_fgINT8(FLGR_Data1D *datdest, FLGR_Data1D *datfft) {
  flgr_no_define_type_function(FLGR_INT8);
}
void flgr1d_fft_backward_cplx_fgINT16(FLGR_Data1D *datdest, FLGR_Data1D *datfft) {
  flgr_no_define_type_function(FLGR_INT16);
}
void flgr1d_fft_backward_cplx_fgINT32(FLGR_Data1D *datdest, FLGR_Data1D *datfft) {
  flgr_no_define_type_function(FLGR_INT32);
}
void flgr1d_fft_backward_cplx_fgFLOAT32(FLGR_Data1D *datdest, FLGR_Data1D *datfft) {
  FLGR_BACKWARD_COMPLEX_FFT_1D(fgFLOAT32,FLGR_FLOAT32,fftwf);
}
void flgr1d_fft_backward_cplx_fgFLOAT64(FLGR_Data1D *datdest, FLGR_Data1D *datfft) {
  FLGR_BACKWARD_COMPLEX_FFT_1D(fgFLOAT64,FLGR_FLOAT64,fftw);
}

////////////////////////////////////////////////////////////////////////////
/*! Compute a Complex backward fft
 *  @param datdest : pointer to FLGR_Data1D
 *  @param datfft  : pointer to FLGR_Data1D
 *  @returns FLGR_RET_OK, ...
 */
//////////////////////////////////////////////////////////////////
FLGR_Ret flgr1d_fft_backward_cplx(FLGR_Data1D *datdest, FLGR_Data1D *datfft) {
  FLGR_Ret ret;

  FLGR_DECLARE_FUNCTION;

  if((datfft==NULL) || (datdest==NULL)) {
    POST_ERROR("Null objects!\n");
    return FLGR_RET_NULL_OBJECT;
  }

  if(datdest->spp!=2)  {
    POST_ERROR("images must be complex (spp=2)\n");
    return FLGR_RET_UNDEFINED_ERROR;
  }

  if(datfft->spp!=2)  {
    POST_ERROR("images must be complex (spp=2)\n");
    return FLGR_RET_UNDEFINED_ERROR;
  }

  if((ret=flgr1d_is_data_same_attributes(datdest,datfft,__FUNCTION__)) != FLGR_RET_OK) return ret;

  FLGR_DISPATCH_PROCEDURE(datdest->type, flgr1d_fft_backward_cplx, datdest, datfft);

  return FLGR_RET_OK;
}



///////////////////////////////////////////////////////////////////
// Macro for different type of flgr1d_fft_backward_cplx_fft_dtype
///////////////////////////////////////////////////////////////////
#define FLGR_MACRO_FFT_FORWARD_REAL_1D(dtype,dtypename,prefix)		\
  FILE *wisdom_file=NULL;						\
  void *bufsrc,*buffft;							\
  prefix##_plan p_forw;							\
  int nx=datsrc->length;						\
  									\
  FLGR_DECLARE_FUNCTION;						\
									\
  bufsrc = flgr_malloc_align(nx*sizeof(dtype),16);			\
  buffft = flgr_malloc_align((nx/2+1)*sizeof(dtype)*2,16);		\
									\
  flgr1d_export_raw_ptr(bufsrc,datsrc);					\
 									\
  wisdom_file=flgr1d_fft_look_for_fftw_wisdom_file(nx,1,dtypename,	\
						   "forward_r2c");	\
  if(wisdom_file!=NULL) prefix##_import_wisdom_from_file(wisdom_file);	\
									\
  p_forw = prefix##_plan_dft_r2c_1d(nx,bufsrc,buffft,FFTW_MEASURE);	\
									\
  if(wisdom_file==NULL) {						\
    wisdom_file=flgr1d_fft_create_fftw_wisdom_file(nx,1,dtypename,	\
						   "forward_r2c");	\
    prefix##_export_wisdom_to_file(wisdom_file);			\
    POST_INFO("fftw wisdom created\n");					\
  }									\
									\
  prefix##_execute(p_forw);						\
									\
  prefix##_destroy_plan(p_forw);					\
									\
  fclose(wisdom_file);							\
									\
  flgr1d_import_raw_##dtype(datfft,buffft);				\
									\
  flgr_free_align(bufsrc);						\
  flgr_free_align(buffft)


void flgr1d_fft_forward_real_fgBIT(FLGR_Data1D *datfft, FLGR_Data1D *datsrc) {
  flgr_no_define_type_function(FLGR_BIT);
}
void flgr1d_fft_forward_real_fgUINT8(FLGR_Data1D *datfft, FLGR_Data1D *datsrc) {
  flgr_no_define_type_function(FLGR_UINT8);
}
void flgr1d_fft_forward_real_fgUINT16(FLGR_Data1D *datfft, FLGR_Data1D *datsrc) {
  flgr_no_define_type_function(FLGR_UINT16);
}
void flgr1d_fft_forward_real_fgUINT32(FLGR_Data1D *datfft, FLGR_Data1D *datsrc) {
  flgr_no_define_type_function(FLGR_UINT32);
}
void flgr1d_fft_forward_real_fgINT8(FLGR_Data1D *datfft, FLGR_Data1D *datsrc) {
  flgr_no_define_type_function(FLGR_INT8);
}
void flgr1d_fft_forward_real_fgINT16(FLGR_Data1D *datfft, FLGR_Data1D *datsrc) {
  flgr_no_define_type_function(FLGR_INT16);
}
void flgr1d_fft_forward_real_fgINT32(FLGR_Data1D *datfft, FLGR_Data1D *datsrc) {
  flgr_no_define_type_function(FLGR_INT32);
}
void flgr1d_fft_forward_real_fgFLOAT32(FLGR_Data1D *datfft, FLGR_Data1D *datsrc) {
  FLGR_MACRO_FFT_FORWARD_REAL_1D(fgFLOAT32,FLGR_FLOAT32,fftwf);
}
void flgr1d_fft_forward_real_fgFLOAT64(FLGR_Data1D *datfft, FLGR_Data1D *datsrc) {
  FLGR_MACRO_FFT_FORWARD_REAL_1D(fgFLOAT64,FLGR_FLOAT64,fftw);
}

////////////////////////////////////////////////////////////////////////////
/*! Compute a Real fft
 *  @param datfft : pointer to FLGR_Data1D
 *  @param datsrc : pointer to FLGR_Data1D
 *  @returns FLGR_RET_OK, ...
 */
//////////////////////////////////////////////////////////////////
FLGR_Ret flgr1d_fft_forward_real(FLGR_Data1D *datfft, FLGR_Data1D *datsrc) {
  FLGR_Ret ret;
  int nx;

  FLGR_DECLARE_FUNCTION;

  if((datsrc==NULL) || (datfft==NULL)){
    POST_ERROR("Null objects!\n");
    return FLGR_RET_NULL_OBJECT;
  }

  nx = datsrc->length;

  if(datsrc->spp!=1)  {
    POST_ERROR("Multi Sample Per Pixel not supported\n");
    return FLGR_RET_NULL_OBJECT;
  }

  if(datfft->spp!=2)  {
    POST_ERROR("data containing fft must be declared as complex (spp=2)\n");
    return FLGR_RET_NULL_OBJECT;
  }

  if((ret=flgr1d_is_data_same_type(datsrc,datfft)) != FLGR_RET_OK) return ret;

  if((nx/2+1) != datfft->length) {
    POST_ERROR("data containing fft must have a size equal to %d\n",nx/2+1);
    return FLGR_RET_SIZE_ERROR;
  }

  FLGR_DISPATCH_PROCEDURE(datsrc->type, flgr1d_fft_forward_real, datfft,datsrc);
}


///////////////////////////////////////////////////////////////////
// Macro for different type of flgr1d_fft_backward_cplx_fft_dtype
///////////////////////////////////////////////////////////////////
#define FLGR_MACRO_FFT_BACKWARD_REAL_1D(dtype,dtypename,prefix)		\
  FILE *wisdom_file=NULL;						\
  void *bufdest, *buffft;						\
  prefix##_plan p_back;							\
  int nx=datdest->length;						\
									\
  FLGR_DECLARE_FUNCTION;						\
									\
  bufdest = flgr_malloc_align(nx*sizeof(dtype),16);			\
  buffft = flgr_malloc_align((nx/2+1)*sizeof(dtype)*2,16);		\
									\
  flgr1d_export_raw_##dtype(buffft,datfft);				\
									\
  wisdom_file=flgr1d_fft_look_for_fftw_wisdom_file(nx,0,dtypename,	\
						   "backward_c2r");	\
  if(wisdom_file!=NULL) prefix##_import_wisdom_from_file(wisdom_file);	\
									\
  p_back = prefix##_plan_dft_c2r_1d(nx,buffft,bufdest,FFTW_MEASURE);	\
									\
  if(wisdom_file==NULL) {						\
    wisdom_file=flgr1d_fft_create_fftw_wisdom_file(nx,0,dtypename,	\
						   "backward_c2r");	\
    prefix##_export_wisdom_to_file(wisdom_file);			\
    POST_INFO("fftw wisdom created\n");					\
  }									\
									\
  prefix##_execute(p_back);						\
									\
  prefix##_destroy_plan(p_back);					\
									\
  fclose(wisdom_file);							\
									\
  flgr1d_import_raw_##dtype(datdest,bufdest);				\
									\
  flgr_free_align(bufdest);						\
  flgr_free_align(buffft)


void flgr1d_fft_backward_real_fgBIT(FLGR_Data1D *datdest, FLGR_Data1D *datfft) {
  flgr_no_define_type_function(FLGR_BIT);
}
void flgr1d_fft_backward_real_fgUINT8(FLGR_Data1D *datdest, FLGR_Data1D *datfft) {
  flgr_no_define_type_function(FLGR_UINT8);
}
void flgr1d_fft_backward_real_fgUINT16(FLGR_Data1D *datdest, FLGR_Data1D *datfft) {
  flgr_no_define_type_function(FLGR_UINT16);
}
void flgr1d_fft_backward_real_fgUINT32(FLGR_Data1D *datdest, FLGR_Data1D *datfft) {
  flgr_no_define_type_function(FLGR_UINT32);
}
void flgr1d_fft_backward_real_fgINT8(FLGR_Data1D *datdest, FLGR_Data1D *datfft) {
  flgr_no_define_type_function(FLGR_INT8);
}
void flgr1d_fft_backward_real_fgINT16(FLGR_Data1D *datdest, FLGR_Data1D *datfft) {
  flgr_no_define_type_function(FLGR_INT16);
}
void flgr1d_fft_backward_real_fgINT32(FLGR_Data1D *datdest, FLGR_Data1D *datfft) {
  flgr_no_define_type_function(FLGR_INT32);
}
void flgr1d_fft_backward_real_fgFLOAT32(FLGR_Data1D *datdest, FLGR_Data1D *datfft) {
  FLGR_MACRO_FFT_BACKWARD_REAL_1D(fgFLOAT32,FLGR_FLOAT32,fftwf);
}
void flgr1d_fft_backward_real_fgFLOAT64(FLGR_Data1D *datdest, FLGR_Data1D *datfft) {
  FLGR_MACRO_FFT_BACKWARD_REAL_1D(fgFLOAT64,FLGR_FLOAT64,fftw);
}


////////////////////////////////////////////////////////////////////////////
/*! Compute a Real backward fft
 *  @param datdest : pointer to FLGR_Data1D
 *  @param datfft  : pointer to FLGR_Data1D
 *  @returns FLGR_RET_OK, ...
 */
//////////////////////////////////////////////////////////////////
FLGR_Ret flgr1d_fft_backward_real(FLGR_Data1D *datdest, FLGR_Data1D *datfft) {
  FLGR_Ret ret;
  int nx;

  FLGR_DECLARE_FUNCTION;

  if((datdest==NULL) || (datfft==NULL)){
    POST_ERROR("Null objects!\n");
    return FLGR_RET_NULL_OBJECT;
  }

  nx = datdest->length;

  if(datdest->spp!=1)  {
    POST_ERROR("Multi Sample Per Pixel not supported\n");
    return FLGR_RET_UNDEFINED_ERROR;
  }

  if(datfft->spp!=2)  {
    POST_ERROR("data containing fft must be declared as complex (spp=2)\n");
    return FLGR_RET_NULL_OBJECT;
  }

  if((ret=flgr1d_is_data_same_type(datdest,datfft)) != FLGR_RET_OK) return ret;

  if((nx/2+1) != datfft->length) {
    POST_ERROR("data containing fft must have a size equal to %d\n",nx/2+1);
    return FLGR_RET_SIZE_ERROR;
  }

  FLGR_DISPATCH_PROCEDURE(datdest->type, flgr1d_fft_backward_real, datdest,datfft);

  return FLGR_RET_OK;
}











/**************************************** 2D Functions ***************************/





























///////////////////////////////////////////////////////////////////
// Macro for different type of flgr2d_fft_forward_cplx
///////////////////////////////////////////////////////////////////
#define FLGR_FORWARD_COMPLEX_FFT_2D(dtype,dtypename,prefix)		\
  FILE *wisdom_file=NULL;						\
  void *bufsrc,*buffft;							\
  prefix##_plan p_forw;							\
  int nx=dat->size_x;							\
  int ny=dat->size_y;							\
  									\
  FLGR_DECLARE_FUNCTION;						\
									\
  bufsrc = fftw_malloc(nx*ny*sizeof(prefix##_complex));			\
  buffft = fftw_malloc(nx*ny*sizeof(prefix##_complex));			\
									\
  flgr2d_export_raw_##dtype(bufsrc,dat);				\
									\
  wisdom_file=flgr2d_fft_look_for_fftw_wisdom_file(ny,nx,1,dtypename,	\
						   "forward_c2c");	\
  if(wisdom_file!=NULL) prefix##_import_wisdom_from_file(wisdom_file);	\
									\
  p_forw = prefix##_plan_dft_2d(ny, nx, bufsrc, buffft,			\
				FFTW_FORWARD,FFTW_MEASURE);		\
  									\
  if(wisdom_file==NULL) {						\
    wisdom_file=flgr2d_fft_create_fftw_wisdom_file(ny,nx,1,dtypename,	\
						   "forward_c2c");	\
    prefix##_export_wisdom_to_file(wisdom_file);			\
    POST_INFO("fftw wisdom created\n");					\
  }									\
									\
  prefix##_execute(p_forw);						\
									\
  prefix##_destroy_plan(p_forw);					\
									\
  fclose(wisdom_file);							\
  									\
  flgr2d_import_raw_##dtype(datfft,buffft);				\
									\
  fftw_free(bufsrc);							\
  fftw_free(buffft)


void flgr2d_fft_forward_cplx_fgBIT(FLGR_Data2D *datfft, FLGR_Data2D *dat) {
  flgr_no_define_type_function(FLGR_BIT);
}
void flgr2d_fft_forward_cplx_fgUINT8(FLGR_Data2D *datfft, FLGR_Data2D *dat) {
  flgr_no_define_type_function(FLGR_UINT8);
}
void flgr2d_fft_forward_cplx_fgUINT16(FLGR_Data2D *datfft, FLGR_Data2D *dat) {
  flgr_no_define_type_function(FLGR_UINT16);
}
void flgr2d_fft_forward_cplx_fgUINT32(FLGR_Data2D *datfft, FLGR_Data2D *dat) {
  flgr_no_define_type_function(FLGR_UINT32);
}
void flgr2d_fft_forward_cplx_fgINT8(FLGR_Data2D *datfft, FLGR_Data2D *dat) {
  flgr_no_define_type_function(FLGR_INT8);
}
void flgr2d_fft_forward_cplx_fgINT16(FLGR_Data2D *datfft, FLGR_Data2D *dat) {
  flgr_no_define_type_function(FLGR_INT16);
}
void flgr2d_fft_forward_cplx_fgINT32(FLGR_Data2D *datfft, FLGR_Data2D *dat) {
  flgr_no_define_type_function(FLGR_INT32);
}
void flgr2d_fft_forward_cplx_fgFLOAT32(FLGR_Data2D *datfft, FLGR_Data2D *dat) {
  FLGR_FORWARD_COMPLEX_FFT_2D(fgFLOAT32,FLGR_FLOAT32,fftwf);
}
void flgr2d_fft_forward_cplx_fgFLOAT64(FLGR_Data2D *datfft, FLGR_Data2D *dat) {
  FLGR_FORWARD_COMPLEX_FFT_2D(fgFLOAT64,FLGR_FLOAT64,fftw);
}

////////////////////////////////////////////////////////////////////////////
/*! Compute a Complex fft
 *  @param datfft : pointer to FLGR_Data2D complex (spp=2)
 *  @param datsrc : pointer to FLGR_Data2D complex (spp=2)
 *  @returns FLGR_RET_OK, ...
 */
//////////////////////////////////////////////////////////////////
FLGR_Ret flgr2d_fft_forward_cplx(FLGR_Data2D *datfft, FLGR_Data2D *datsrc) {
  FLGR_Ret ret;

  FLGR_DECLARE_FUNCTION;

  if((datsrc == NULL) || (datfft == NULL)) {
    POST_ERROR("Null objects!\n");
    return FLGR_RET_NULL_OBJECT;
  }

  if(datsrc->spp!=2)  {
    POST_ERROR("Only complex images are supported\n");
    return FLGR_RET_NULL_OBJECT;
  }

  if((ret=flgr2d_is_data_same_attributes(datsrc,datfft,__FUNCTION__)) != FLGR_RET_OK) return ret;

  FLGR_DISPATCH_PROCEDURE(datsrc->type, flgr2d_fft_forward_cplx, datfft,datsrc);
}




///////////////////////////////////////////////////////////////////
// Macro for different type of flgr2d_fft_backward_cplx
///////////////////////////////////////////////////////////////////
#define FLGR_BACKWARD_COMPLEX_FFT_2D(dtype,dtypename,prefix)		\
  FILE *wisdom_file=NULL;						\
  void *bufdest, *buffft;						\
  prefix##_plan p_back;							\
  int nx=datfft->size_x;						\
  int ny=datfft->size_y;						\
									\
  FLGR_DECLARE_FUNCTION;						\
									\
  bufdest = fftw_malloc(nx*ny*sizeof(prefix##_complex));		\
  buffft = fftw_malloc(nx*ny*sizeof(prefix##_complex));			\
									\
  flgr2d_export_raw_##dtype(buffft,datfft);				\
									\
  wisdom_file=flgr2d_fft_look_for_fftw_wisdom_file(ny,nx,0,dtypename,	\
						   "backward_c2c");	\
  if(wisdom_file!=NULL) prefix##_import_wisdom_from_file(wisdom_file);	\
									\
  p_back = prefix##_plan_dft_2d(ny, nx, buffft, bufdest,		\
				FFTW_BACKWARD,FFTW_MEASURE);		\
  									\
  if(wisdom_file==NULL) {						\
    wisdom_file=flgr2d_fft_create_fftw_wisdom_file(ny,nx,0,dtypename,	\
						   "backward_c2c");	\
    prefix##_export_wisdom_to_file(wisdom_file);			\
    POST_INFO("fftw wisdom created\n");					\
  }									\
  									\
  prefix##_execute(p_back);						\
  									\
  prefix##_destroy_plan(p_back);					\
  									\
  fclose(wisdom_file);							\
									\
  flgr2d_import_raw_##dtype(datdest,bufdest);				\
									\
  fftw_free(bufdest);							\
  fftw_free(buffft);							\
									\
  return


void flgr2d_fft_backward_cplx_fgBIT(FLGR_Data2D *datdest, FLGR_Data2D *datfft) {
  flgr_no_define_type_function(FLGR_BIT);
}
void flgr2d_fft_backward_cplx_fgUINT8(FLGR_Data2D *datdest, FLGR_Data2D *datfft) {
  flgr_no_define_type_function(FLGR_UINT8);
}
void flgr2d_fft_backward_cplx_fgUINT16(FLGR_Data2D *datdest, FLGR_Data2D *datfft) {
  flgr_no_define_type_function(FLGR_UINT16);
}
void flgr2d_fft_backward_cplx_fgUINT32(FLGR_Data2D *datdest, FLGR_Data2D *datfft) {
  flgr_no_define_type_function(FLGR_UINT32);
}
void flgr2d_fft_backward_cplx_fgINT8(FLGR_Data2D *datdest, FLGR_Data2D *datfft) {
  flgr_no_define_type_function(FLGR_INT8);
}
void flgr2d_fft_backward_cplx_fgINT16(FLGR_Data2D *datdest, FLGR_Data2D *datfft) {
  flgr_no_define_type_function(FLGR_INT16);
}
void flgr2d_fft_backward_cplx_fgINT32(FLGR_Data2D *datdest, FLGR_Data2D *datfft) {
  flgr_no_define_type_function(FLGR_INT32);
}
void flgr2d_fft_backward_cplx_fgFLOAT32(FLGR_Data2D *datdest, FLGR_Data2D *datfft) {
  FLGR_BACKWARD_COMPLEX_FFT_2D(fgFLOAT32,FLGR_FLOAT32,fftwf);
}
void flgr2d_fft_backward_cplx_fgFLOAT64(FLGR_Data2D *datdest, FLGR_Data2D *datfft) {
  FLGR_BACKWARD_COMPLEX_FFT_2D(fgFLOAT64,FLGR_FLOAT64,fftw);
}

////////////////////////////////////////////////////////////////////////////
/*! Compute a Complex backward fft
 *  @param datdest : pointer to FLGR_Data2D
 *  @param datfft  : pointer to FLGR_Data2D
 *  @returns FLGR_RET_OK, ...
 */
//////////////////////////////////////////////////////////////////
FLGR_Ret flgr2d_fft_backward_cplx(FLGR_Data2D *datdest, FLGR_Data2D *datfft) {
  FLGR_Ret ret;

  FLGR_DECLARE_FUNCTION;

  if((datfft==NULL) || (datdest==NULL)) {
    POST_ERROR("Null objects!\n");
    return FLGR_RET_NULL_OBJECT;
  }

  if(datdest->spp!=2)  {
    POST_ERROR("images must be complex (spp=2)\n");
    return FLGR_RET_UNDEFINED_ERROR;
  }

  if(datfft->spp!=2)  {
    POST_ERROR("images must be complex (spp=2)\n");
    return FLGR_RET_UNDEFINED_ERROR;
  }

  if((ret=flgr2d_is_data_same_attributes(datdest,datfft,__FUNCTION__)) != FLGR_RET_OK) return ret;

  FLGR_DISPATCH_PROCEDURE(datdest->type, flgr2d_fft_backward_cplx, datdest, datfft);

  return FLGR_RET_OK;
}



///////////////////////////////////////////////////////////////////
// Macro for different type of flgr2d_fft_backward_cplx_fft_dtype
///////////////////////////////////////////////////////////////////
#define FLGR_MACRO_FFT_FORWARD_REAL_2D(dtype,dtypename,prefix)		\
  FILE *wisdom_file=NULL;						\
  void *bufsrc,*buffft;							\
  prefix##_plan p_forw;							\
  int nx=datsrc->size_x;						\
  int ny=datsrc->size_y;						\
									\
  FLGR_DECLARE_FUNCTION;						\
									\
  bufsrc = flgr_malloc_align(nx*ny*sizeof(dtype),16);			\
  buffft = flgr_malloc_align((nx/2+1)*ny*sizeof(dtype)*2,16);		\
									\
  flgr2d_export_raw_ptr(bufsrc,datsrc);					\
 									\
  wisdom_file=flgr2d_fft_look_for_fftw_wisdom_file(ny,nx,1,dtypename,	\
						   "forward_r2c");	\
  if(wisdom_file!=NULL) prefix##_import_wisdom_from_file(wisdom_file);	\
									\
  p_forw = prefix##_plan_dft_r2c_2d(ny, nx,bufsrc,buffft,FFTW_MEASURE);	\
									\
  if(wisdom_file==NULL) {						\
    wisdom_file=flgr2d_fft_create_fftw_wisdom_file(ny,nx,1,dtypename,	\
						   "forward_r2c");	\
    prefix##_export_wisdom_to_file(wisdom_file);			\
    POST_INFO("fftw wisdom created\n");					\
  }									\
									\
  prefix##_execute(p_forw);						\
									\
  prefix##_destroy_plan(p_forw);					\
									\
  fclose(wisdom_file);							\
									\
  flgr2d_import_raw_##dtype(datfft,buffft);				\
									\
  fftw_free(bufsrc);							\
  fftw_free(buffft)


void flgr2d_fft_forward_real_fgBIT(FLGR_Data2D *datfft, FLGR_Data2D *datsrc) {
  flgr_no_define_type_function(FLGR_BIT);
}
void flgr2d_fft_forward_real_fgUINT8(FLGR_Data2D *datfft, FLGR_Data2D *datsrc) {
  flgr_no_define_type_function(FLGR_UINT8);
}
void flgr2d_fft_forward_real_fgUINT16(FLGR_Data2D *datfft, FLGR_Data2D *datsrc) {
  flgr_no_define_type_function(FLGR_UINT16);
}
void flgr2d_fft_forward_real_fgUINT32(FLGR_Data2D *datfft, FLGR_Data2D *datsrc) {
  flgr_no_define_type_function(FLGR_UINT32);
}
void flgr2d_fft_forward_real_fgINT8(FLGR_Data2D *datfft, FLGR_Data2D *datsrc) {
  flgr_no_define_type_function(FLGR_INT8);
}
void flgr2d_fft_forward_real_fgINT16(FLGR_Data2D *datfft, FLGR_Data2D *datsrc) {
  flgr_no_define_type_function(FLGR_INT16);
}
void flgr2d_fft_forward_real_fgINT32(FLGR_Data2D *datfft, FLGR_Data2D *datsrc) {
  flgr_no_define_type_function(FLGR_INT32);
}
void flgr2d_fft_forward_real_fgFLOAT32(FLGR_Data2D *datfft, FLGR_Data2D *datsrc) {
  FLGR_MACRO_FFT_FORWARD_REAL_2D(fgFLOAT32,FLGR_FLOAT32,fftwf);
}
void flgr2d_fft_forward_real_fgFLOAT64(FLGR_Data2D *datfft, FLGR_Data2D *datsrc) {
  FLGR_MACRO_FFT_FORWARD_REAL_2D(fgFLOAT64,FLGR_FLOAT64,fftw);
}

////////////////////////////////////////////////////////////////////////////
/*! Compute a Real fft
 *  @param datsrc : pointer to FLGR_Data2D
 *  @param datfft : pointer to FLGR_Data2D
 *  @returns FLGR_RET_OK, ...
 */
//////////////////////////////////////////////////////////////////
FLGR_Ret flgr2d_fft_forward_real(FLGR_Data2D *datfft, FLGR_Data2D *datsrc) {
  FLGR_Ret ret;
  int nx;

  FLGR_DECLARE_FUNCTION;

  if((datsrc==NULL) || (datfft==NULL)){
    POST_ERROR("Null objects!\n");
    return FLGR_RET_NULL_OBJECT;
  }

  nx = datsrc->size_x;

  if(datsrc->spp!=1)  {
    POST_ERROR("Multi Sample Per Pixel not supported\n");
    return FLGR_RET_NULL_OBJECT;
  }

  if(datfft->spp!=2)  {
    POST_ERROR("data containing fft must be complex (spp=2)\n");
    return FLGR_RET_NULL_OBJECT;
  }

  if((ret=flgr2d_is_data_same_type(datsrc,datfft)) != FLGR_RET_OK) return ret;

  if((nx/2+1) != datfft->size_x) {
    POST_ERROR("data containing fft must have a width equal to %d\n",nx/2+1);
    return FLGR_RET_SIZE_ERROR;
  }

 
  FLGR_DISPATCH_PROCEDURE(datsrc->type, flgr2d_fft_forward_real, datfft,datsrc);
}


///////////////////////////////////////////////////////////////////
// Macro for different type of flgr2d_fft_backward_cplx_fft_dtype
///////////////////////////////////////////////////////////////////
#define FLGR_MACRO_FFT_BACKWARD_REAL_2D(dtype,dtypename,prefix)		\
  FILE *wisdom_file=NULL;						\
  void *bufdest, *buffft;						\
  prefix##_plan p_back;							\
  int nx=datdest->size_x;						\
  int ny=datdest->size_y;						\
									\
  FLGR_DECLARE_FUNCTION;						\
									\
  bufdest = flgr_malloc_align(nx*ny*sizeof(dtype),16);			\
  buffft = flgr_malloc_align((nx/2+1)*ny*sizeof(dtype)*2,16);		\
									\
  flgr2d_export_raw_##dtype(buffft,datfft);				\
									\
  wisdom_file=flgr2d_fft_look_for_fftw_wisdom_file(ny,nx,0,dtypename,	\
						   "backward_c2r");	\
  if(wisdom_file!=NULL) prefix##_import_wisdom_from_file(wisdom_file);	\
									\
  p_back = prefix##_plan_dft_c2r_2d(ny,nx,buffft,bufdest,FFTW_MEASURE);	\
									\
  if(wisdom_file==NULL) {						\
    wisdom_file=flgr2d_fft_create_fftw_wisdom_file(ny,nx,0,dtypename,	\
						   "backward_c2r");	\
    prefix##_export_wisdom_to_file(wisdom_file);			\
    POST_INFO("fftw wisdom created\n");					\
  }									\
									\
  prefix##_execute(p_back);						\
									\
  prefix##_destroy_plan(p_back);					\
									\
  fclose(wisdom_file);							\
									\
  flgr2d_import_raw_##dtype(datdest,bufdest);				\
									\
  flgr_free_align(bufdest);						\
  flgr_free_align(buffft)


void flgr2d_fft_backward_real_fgBIT(FLGR_Data2D *datdest, FLGR_Data2D *datfft) {
  flgr_no_define_type_function(FLGR_BIT);
}
void flgr2d_fft_backward_real_fgUINT8(FLGR_Data2D *datdest, FLGR_Data2D *datfft) {
  flgr_no_define_type_function(FLGR_UINT8);
}
void flgr2d_fft_backward_real_fgUINT16(FLGR_Data2D *datdest, FLGR_Data2D *datfft) {
  flgr_no_define_type_function(FLGR_UINT16);
}
void flgr2d_fft_backward_real_fgUINT32(FLGR_Data2D *datdest, FLGR_Data2D *datfft) {
  flgr_no_define_type_function(FLGR_UINT32);
}
void flgr2d_fft_backward_real_fgINT8(FLGR_Data2D *datdest, FLGR_Data2D *datfft) {
  flgr_no_define_type_function(FLGR_INT8);
}
void flgr2d_fft_backward_real_fgINT16(FLGR_Data2D *datdest, FLGR_Data2D *datfft) {
  flgr_no_define_type_function(FLGR_INT16);
}
void flgr2d_fft_backward_real_fgINT32(FLGR_Data2D *datdest, FLGR_Data2D *datfft) {
  flgr_no_define_type_function(FLGR_INT32);
}
void flgr2d_fft_backward_real_fgFLOAT32(FLGR_Data2D *datdest, FLGR_Data2D *datfft) {
  FLGR_MACRO_FFT_BACKWARD_REAL_2D(fgFLOAT32,FLGR_FLOAT32,fftwf);
}
void flgr2d_fft_backward_real_fgFLOAT64(FLGR_Data2D *datdest, FLGR_Data2D *datfft) {
  FLGR_MACRO_FFT_BACKWARD_REAL_2D(fgFLOAT64,FLGR_FLOAT64,fftw);
}


////////////////////////////////////////////////////////////////////////////
/*! Compute a Real backward fft
 *  @param datdest : pointer to FLGR_Data2D
 *  @param datfft : pointer to FLGR_Data2D
 *  @returns FLGR_RET_OK, ...
 */
//////////////////////////////////////////////////////////////////
FLGR_Ret flgr2d_fft_backward_real(FLGR_Data2D *datdest, FLGR_Data2D *datfft) {
  FLGR_Ret ret;
  int nx;

  FLGR_DECLARE_FUNCTION;

  if((datdest==NULL) || (datfft==NULL)){
    POST_ERROR("Null objects!\n");
    return FLGR_RET_NULL_OBJECT;
  }

  nx = datdest->size_x;

  if(datdest->spp!=1)  {
    POST_ERROR("Multi Sample Per Pixel not supported\n");
    return FLGR_RET_UNDEFINED_ERROR;
  }


  if(datfft->spp!=2)  {
    POST_ERROR("data containing fft must be complex (spp=2)\n");
    return FLGR_RET_NULL_OBJECT;
  }

  if((ret=flgr2d_is_data_same_type(datdest,datfft)) != FLGR_RET_OK) return ret;

  if((nx/2+1) != datfft->size_x) {
    POST_ERROR("data containing fft must have a width equal to %d\n",nx/2+1);
    return FLGR_RET_SIZE_ERROR;
  }

  FLGR_DISPATCH_PROCEDURE(datdest->type, flgr2d_fft_backward_real, datdest,datfft);

  return FLGR_RET_OK;
}








///////////////////////////////////////////////////////////////////
// Macro for different type of switchCorner_dtype
///////////////////////////////////////////////////////////////////
#define FLGR_MACRO_SWITCH_CORNER_2D(dtype)	\
  dtype *pdestA,*pdestB,*pdestC,*pdestD;	\
  dtype *psrcA,*psrcB,*psrcC,*psrcD;		\
  int i;					\
  int spp = datdest->spp;			\
  int h2 = datdest->size_y>>1;			\
  int w2 = (datdest->size_x>>1)*spp;		\
 						\
  FLGR_DECLARE_FUNCTION;			\
						\
  for(i=0 ; i < h2 ; i++) {			\
    psrcA = (dtype *) datsrc->array[i];		\
    psrcB = psrcA+w2;				\
    psrcC = (dtype *) datsrc->array[i+h2];	\
    psrcD = psrcC+w2;				\
						\
    pdestB = (dtype *) datdest->array[i+h2];	\
    pdestA = pdestB+w2;				\
    pdestD = (dtype *) datdest->array[i];	\
    pdestC = pdestD+w2;				\
						\
    memcpy(pdestA,psrcA,sizeof(dtype)*w2);	\
    memcpy(pdestB,psrcB,sizeof(dtype)*w2);	\
    memcpy(pdestC,psrcC,sizeof(dtype)*w2);	\
    memcpy(pdestD,psrcD,sizeof(dtype)*w2);	\
						\
  }						\
						\
  return



void flgr2d_fft_switch_corner_fgBIT(FLGR_Data2D *datdest, FLGR_Data2D *datsrc) {
  flgr_no_define_type_function(FLGR_BIT);
}
void flgr2d_fft_switch_corner_fgUINT8(FLGR_Data2D *datdest, FLGR_Data2D *datsrc) {
  FLGR_MACRO_SWITCH_CORNER_2D(fgUINT8);
}
void flgr2d_fft_switch_corner_fgUINT16(FLGR_Data2D *datdest, FLGR_Data2D *datsrc) {
  FLGR_MACRO_SWITCH_CORNER_2D(fgUINT16);
}
void flgr2d_fft_switch_corner_fgUINT32(FLGR_Data2D *datdest, FLGR_Data2D *datsrc) {
  FLGR_MACRO_SWITCH_CORNER_2D(fgUINT32);
}
void flgr2d_fft_switch_corner_fgINT8(FLGR_Data2D *datdest, FLGR_Data2D *datsrc) {
  FLGR_MACRO_SWITCH_CORNER_2D(fgINT8);
}
void flgr2d_fft_switch_corner_fgINT16(FLGR_Data2D *datdest, FLGR_Data2D *datsrc) {
  FLGR_MACRO_SWITCH_CORNER_2D(fgINT16);
}
void flgr2d_fft_switch_corner_fgINT32(FLGR_Data2D *datdest, FLGR_Data2D *datsrc) {
  FLGR_MACRO_SWITCH_CORNER_2D(fgINT32);
}
void flgr2d_fft_switch_corner_fgFLOAT32(FLGR_Data2D *datdest, FLGR_Data2D *datsrc) {
  FLGR_MACRO_SWITCH_CORNER_2D(fgFLOAT32);
}
void flgr2d_fft_switch_corner_fgFLOAT64(FLGR_Data2D *datdest, FLGR_Data2D *datsrc) {
  FLGR_MACRO_SWITCH_CORNER_2D(fgFLOAT64);
}
////////////////////////////////////////////////////////////////////////////
/*! Sort the 4 corners after a fft
 *  @param datsrc : pointer to FLGR_Data2D
 *  @param datdest : pointer to FLGR_Data2D
 *  @returns FLGR_RET_OK, ...
 */
//////////////////////////////////////////////////////////////////
FLGR_Ret flgr2d_fft_switch_corner(FLGR_Data2D *datdest, FLGR_Data2D *datsrc) {
  FLGR_Ret ret;

  FLGR_DECLARE_FUNCTION;
  
  if((datdest==NULL) || (datsrc==NULL)){
    POST_ERROR("Null objects!\n");
    return FLGR_RET_NULL_OBJECT;
  }

  if((ret=flgr2d_is_data_same_attributes(datdest,datsrc,__FUNCTION__)) != FLGR_RET_OK) return ret;

  FLGR_DISPATCH_PROCEDURE(datdest->type, flgr2d_fft_switch_corner, datdest,datsrc);

}



//! @}
